home *** CD-ROM | disk | FTP | other *** search
/ CDUTIL 13 / CDUTIL #13 Julio 1995.iso / windows / acadcom / acrx / sample / gravity.cc < prev    next >
Encoding:
Text File  |  1995-02-08  |  67.1 KB  |  1,870 lines

  1. /* Next available MSG number is 141 */
  2. /*
  3.  
  4.    GRAVITY.CC
  5.  
  6.    Copyright (C) 1994 by Autodesk, Inc.
  7.  
  8.    Permission to use, copy, modify, and distribute this software in 
  9.    object code form for any purpose and without fee is hereby granted, 
  10.    provided that the above copyright notice appears in all copies and 
  11.    that both that copyright notice and the limited warranty and 
  12.    restricted rights notice below appear in all supporting 
  13.    documentation.
  14.  
  15.    AUTODESK PROVIDES THIS PROGRAM "AS IS" AND WITH ALL FAULTS.  
  16.    AUTODESK SPECIFICALLY DISCLAIMS ANY IMPLIED WARRANTY OF 
  17.    MERCHANTABILITY OR FITNESS FOR A PARTICULAR USE.  AUTODESK, INC.
  18.    DOES NOT WARRANT THAT THE OPERATION OF THE PROGRAM WILL BE 
  19.    UNINTERRUPTED OR ERROR FREE.
  20.  
  21.    Use, duplication, or disclosure by the U.S. Government is subject to 
  22.    restrictions set forth in FAR 52.227-19 (Commercial Computer 
  23.    Software - Restricted Rights) and DFAR 252.227-7013(c)(1)(ii) 
  24.    (Rights in Technical Data and Computer Software), as applicable.
  25.     
  26.    .
  27.  
  28.  
  29.         N-Body Gravitational Interaction Simulator for AutoCAD
  30.  
  31.         A sample ADS application.
  32.  
  33.         Designed and implemented by John Walker in August of 1989.
  34.  
  35.            "Ubi materia, ibi geometria."
  36.                       -- Johannes Kepler
  37.  
  38.            "The orbit of  any  one  planet  depends  on  the
  39.            combined  motion  of  all  the  planets,  not  to
  40.            mention the action of all of these on each other.
  41.            But  to  consider simultaneously all these causes
  42.            of motion and to define these  motions  by  exact
  43.            laws   allowing   of   conventional   calculation
  44.            exceeds, unless I am mistaken, the force  of  the
  45.            entire human intellect."
  46.                       -- Sir Isaac Newton, Principia
  47.  
  48.  
  49.         The  physics  underlying  this simulation are explained in the
  50.         the  chapter  "A  Cosmic  Ballet",  pp.   229-238  of  A.   K.
  51.         Dewdney,  "The  Armchair Universe", W.  H.  Freeman: New York,
  52.         1988.
  53.  
  54.         The following commands are implemented in this file:
  55.  
  56.         DEMO        Creates  one  of  a series of standard demo models
  57.                     when executed in  a  new,  empty  drawing.   These
  58.                     models  may be modified with the MASS and MASSEDIT
  59.                     commands just like models entered manually.
  60.  
  61.         FRAME       Asks you to pick a mass entity.   Motion  is  then
  62.                     displayed in that mass's reference frame (in other
  63.                     words, that mass  appears  stationary  and  others
  64.                     move  around  it).  The default reference frame is
  65.                     "Inertial", an unaccelerated frame  at  rest  with
  66.                     respect  to  the  distant  stars.   Specifying  no
  67.                     entity to the  FRAME  command  re-establishes  the
  68.                     inertial frame.
  69.  
  70.         MASS        Creates  a  new  mass.   You're invited to enter a
  71.                     name for the mass, its position (specified by  any
  72.                     means  of  co-ordinate  specification), a velocity
  73.                     vector in units of  astronomical  units  per  year
  74.                     (which  can  be either entered explicitly or drawn
  75.                     by typing "V", then dragging the endpoint  of  the
  76.                     velocity  vector to the correct position), and its
  77.                     mass in units of Solar masses.
  78.  
  79.         MASSEDIT    Lets you modify the name, velocity, and mass of an
  80.                     existing mass.  Pick a single mass by pointing.  A
  81.                     dialogue is displayed with the properties  of  the
  82.                     mass.   Change  them  as you wish, then pick OK to
  83.                     update the properties of the mass in the database.
  84.                     If you pick Cancel, the mass is not modified.
  85.  
  86.         RESET       When  a  simulation  is  run,  it  halts after the
  87.                     specified number of steps or  simulated  time  has
  88.                     elapsed.    A   subsequent  RUN  command  normally
  89.                     resumes the simulation from the point at which the
  90.                     last  stopped.  RESET erases all motion paths from
  91.                     the screen and sets the  simulation  back  to  the
  92.                     start.    RUN   after   a  RESET  will  begin  the
  93.                     simulation from the initial state.
  94.  
  95.         RUN         Starts the simulation.  You're  asked  to  specify
  96.                     the length of the simulation as either a number of
  97.                     motion steps or by the number of simulated  years.
  98.                     If  you  enter  a positive number, it's taken as a
  99.                     step count.  Negative numbers (which  may  include
  100.                     decimal   fractions)  specify  simulated  time  in
  101.                     years.  The length in time of each simulation step
  102.                     varies  based  upon the velocity and separation of
  103.                     masses, so if objects approach  one  another  very
  104.                     closely  a  very  large  number  of  steps will be
  105.                     required to simulate a given time  interval.   The
  106.                     calculation time per step is essentially constant,
  107.                     so if you're investigating a system  with  unknown
  108.                     behaviour  you'll probably want to RUN for a given
  109.                     number of cycles, but when running a stable system
  110.                     such as the Solar System, you can simulate a given
  111.                     time span.  In  any  case,  you  can  terminate  a
  112.                     simulation with the Control C key.  RUN  can  also
  113.                     be  invoked as a function, (C:RUN <length>), where
  114.                     <length> specifies  the  simulation  length  by  a
  115.                     positive  or  negative  number as described above.
  116.                     When called  as  a  function,  C:RUN  returns  the
  117.                     simulation end time as its result.
  118.  
  119.         SETGRAV     The  Gravity  simulator has several variables that
  120.                     control its operation.  These variables are  saved
  121.                     with the drawing and may be inspected and modified
  122.                     with the  SETGRAV  command.   SETGRAV  displays  a
  123.                     dialogue  in  which  you  may  change  any  of the
  124.                     following variables:
  125.  
  126.                        Output to display only?
  127.                           This is set to Yes or  No  (True/False,  and
  128.                           1/0  are  also  accepted).   If  set  to the
  129.                           default value of Yes, motion paths are drawn
  130.                           using  temporary  vectors within the current
  131.                           viewport.  If the picture is  REDRAWn,  they
  132.                           disappear.   If  set to No, motion paths are
  133.                           added to the drawing as LINE entities.  This
  134.                           causes paths to appear in all viewports, and
  135.                           you can  ZOOM  on  sections  of  a  path  to
  136.                           examine   it   in   greater  detail.   Paths
  137.                           represented as  LINEs  are  saved  with  the
  138.                           drawing.  Generating LINE entities for paths
  139.                           is much slower than just drawing them on the
  140.                           screen,  so  choose  this mode only when you
  141.                           need it.
  142.  
  143.                        Display step number?
  144.                           If this mode is  "Yes"  (the  default),  the
  145.                           simulation  step  number is displayed in the
  146.                           coordinate status  line  as  the  simulation
  147.                           progresses.
  148.  
  149.                        Display time?
  150.                           If "Yes" (the default) the simulated time in
  151.                           years is displayed in the coordinate  status
  152.                           line.
  153.  
  154.                        Step size?
  155.                           This real number specifies the  factor  used
  156.                           to  determine  the  size  of the integration
  157.                           steps used in calculating the motion of  the
  158.                           masses.   The  time  step  is  calculated by
  159.                           dividing  the  distance  between   the   two
  160.                           closest   masses  by  the  highest  relative
  161.                           velocity  of  any  pair  of  masses.    This
  162.                           quantity,  measured  in years, is multiplied
  163.                           by the factor  given  by  this  variable  to
  164.                           obtain  the length of the step.  The default
  165.                           value  of  0.1  works  well  for  reasonably
  166.                           well-behaved  simulations.  If you find that
  167.                           a simulation is taking  too  long,  you  can
  168.                           speed it up by increasing the step size, but
  169.                           be aware that the results you see may not be
  170.                           physically  accurate.   Increasing  the step
  171.                           size magnifies the inaccuracies of  modeling
  172.                           a  continuous process such as gravitation by
  173.                           discrete  steps.   In   general,   you   can
  174.                           increase the step size as long as you obtain
  175.                           the same results.  When the  outcome  begins
  176.                           to vary, you've set the step size too large.
  177.  
  178.                        Minimum step?
  179.                           After  the  step  size  is   calculated   as
  180.                           described  above  it  is  compared  with the
  181.                           minimum step size and, if less, the  minimum
  182.                           is  used.   The  minimum step size is set to
  183.                           0.00001 years (about  5  minutes);  this  is
  184.                           sufficient  to unstick many simulations that
  185.                           involve a close encounter.  Amazingly,  even
  186.                           a  minimum  this  short  can  yield  grossly
  187.                           inaccurate results when solar masses execute
  188.                           hairpin turns about one another
  189.  
  190.         UPDATE      A  simulation  normally  proceeds from the initial
  191.                     positions, velocities, and masses of  the  objects
  192.                     specified  by  the  MASS  command.  The simulation
  193.                     keeps track of these values as it  progresses  but
  194.                     does  not automatically adjust the mass objects in
  195.                     the database.  If you want to  move  the  database
  196.                     objects  to  the  positions at the end of the most
  197.                     recent simulation (and adjust their velocities  to
  198.                     the  corresponding  instantaneous velocities), use
  199.                     the UPDATE command.  If you don't  do  an  UPDATE,
  200.                     the masses will remain in their original positions
  201.                     even if you  save  the  drawing  after  running  a
  202.                     simulation.
  203.  
  204. */
  205.  
  206. #include   <stdio.h>
  207. #include   <string.h>
  208. #ifndef Macintosh
  209. #include   <ctype.h>
  210. #endif
  211. #include <stdlib.h>
  212. #include   <math.h>
  213. #include    <assert.h>
  214. #include "rxdefs.h"
  215. #include   "adslib.h"
  216.  
  217. /*  Standard drawing object names  */
  218.  
  219. /* Utility frozen layer for information */
  220. #define FrozenLayer /*MSG1*/"FROZEN-SOLID"
  221. #define OrbitLayer  /*MSG2*/"ORBITS"  /* Layer for orbital path traces */
  222.  
  223. /*  Data types  */
  224.  
  225. typedef enum {False = 0, True = 1} Boolean;
  226.  
  227. #define HANDLEN  18                   /* String long enough to hold a handle */
  228.  
  229. /* Definitions  to  wrap  around  submission  of  AutoCAD commands to
  230.    prevent their being echoed.  */
  231.  
  232. #define Cmdecho  False                /* Make True for debug command output */
  233.  
  234. #define CommandB()  { struct resbuf rBc, rBb, rBu, rBh; \
  235.         ads_getvar(/*MSG0*/"CMDECHO", &rBc); \
  236.         ads_getvar(/*MSG0*/"BLIPMODE", &rBb); \
  237.         ads_getvar(/*MSG0*/"HIGHLIGHT", &rBh); \
  238.         rBu.restype = RTSHORT; \
  239.         rBu.resval.rint = (int) Cmdecho; \
  240.         ads_setvar(/*MSG0*/"CMDECHO", &rBu); \
  241.         rBu.resval.rint = (int) False; \
  242.         ads_setvar(/*MSG0*/"BLIPMODE", &rBu); \
  243.         ads_setvar(/*MSG0*/"HIGHLIGHT", &rBu)
  244.  
  245. #define CommandE()  ads_setvar(/*MSG0*/"CMDECHO", &rBc); \
  246.                     ads_setvar(/*MSG0*/"BLIPMODE", &rBb); \
  247.                     ads_setvar(/*MSG0*/"HIGHLIGHT", &rBh); }
  248.  
  249. /*  Definitions  that permit you to push and pop system variables with
  250.     minimal complexity.  These don't work  (and  will  cause  horrible
  251.     crashes  if  used)  with  string  variables,  but since all string
  252.     variables are read-only, they cannot be saved and restored in  any
  253.     case.  */
  254.  
  255. #define PushVar(var, cell, newval, newtype) { struct resbuf cell, cNeW; \
  256.         ads_getvar(var, &cell); cNeW.restype = cell.restype;             \
  257.         cNeW.resval.newtype = newval; ads_setvar(var, &cNeW)
  258.  
  259. #define PopVar(var, cell) ads_setvar(var, &cell); }
  260.  
  261. /* Set point variable from three co-ordinates */
  262.  
  263. #define Spoint(pt, x, y, z)  pt[X] = (x);  pt[Y] = (y);  pt[Z] = (z)
  264.  
  265. /* Copy point  */
  266.  
  267. #define Cpoint(d, s)   d[X] = s[X];  d[Y] = s[Y];  d[Z] = s[Z]
  268.  
  269. /* Generation parameters for objects */
  270.  
  271. #define SphereLong  12                /* Longitudinal tabulations on sphere */
  272. #define SphereLat   12                /* Latitudinal tabulations on sphere */
  273.  
  274. /* Utility definition to get an  array's  element  count  (at  compile
  275.    time).   For  example:
  276.  
  277.        int  arr[] = {1,2,3,4,5};
  278.        ...
  279.        printf("%d", ELEMENTS(arr));
  280.  
  281.    would print a five.  ELEMENTS("abc") can also be used to  tell  how
  282.    many  bytes are in a string constant INCLUDING THE TRAILING NULL. */
  283.  
  284. #define ELEMENTS(array) (sizeof(array)/sizeof((array)[0]))
  285.  
  286. /* Utility definitions */
  287.  
  288. #ifndef abs
  289. #define abs(x)      ((x)<0 ? -(x) : (x))
  290. #endif  /* abs */
  291. #ifdef min
  292. #undef  min
  293. #endif
  294. #define min(a,b)    ((a)<(b) ? (a) : (b))
  295. #ifdef max
  296. #undef  max
  297. #endif
  298. #define max(a,b)    ((a)>(b) ? (a) : (b))
  299.  
  300. /*  Many C implementations may lack a cube root function.  Rather than
  301.     count  on  a system cbrt() function, we define our own in terms of
  302.     functions more likely to be available.  */
  303.  
  304. #define cuberoot(x) (exp(log(x) / 3.0))
  305.  
  306. /* All Function Prototypes for gravity.c */
  307.  
  308. int            main    _((int argc, char *argv[]));
  309. static Boolean  funcload _((void));
  310. static char *   alloc   _((unsigned nbytes));
  311. static void     defmassblk _((void));
  312. static struct resbuf *
  313.                 entitem _((struct resbuf *rchain, int gcode));
  314. static void     entinfo _((ads_name en,char *h,ads_point p,ads_real *r,int *c));
  315. static void     partext _((void));
  316. static void     addvec  _((ads_real *ap, ads_real *bp, ads_real *cp));
  317. static void     subvec  _((ads_real *ap, ads_real *bp, ads_real *cp));
  318. static ads_real sqabsv  _((ads_real *ap));
  319. static void     sumvec  _((ads_real *ap,ads_real *bp,ads_real t, ads_real *cp));
  320. static void     varblockdef _((void));
  321. static void     savemodes _((void));
  322. static Boolean  boolvar _((char *varname, char *s));
  323. static void     varset  _((void));
  324. static void     setframe _((void));
  325. static Boolean  initacad _((Boolean reset));
  326. static void     addmass _((char *mn,ads_point pos,ads_point vel,ads_real mass));
  327. static void     mass    _((void));
  328. static void     demo    _((void));
  329. static void     massedit _((void));
  330. static void     reset   _((void));
  331. static void     frame   _((void));
  332. static void     run     _((void));
  333. static void     setgrav _((void));
  334. static void     update  _((void));
  335. #ifdef  HIGHC
  336. static void     abort   _((void));
  337. #endif
  338.  
  339.  
  340. extern "C" 
  341. {                         
  342. AcRx::AppRetCode acrxEntryPoint(AcRx::AppMsgCode msg,void * pkt);
  343. }
  344.  
  345. /*  This  program  works in a somewhat unconventional system of units.
  346.     Length is measured in astronomical units (the mean  distance  from
  347.     the  Earth  to the Sun), mass in units of the mass of the Sun, and
  348.     time in years.  The following definitions derive the value of  the
  349.     gravitational  constant  in this system of units from its handbook
  350.     definition in SI units.  */
  351.  
  352. #define G_SI        6.6732e-11        /* (Newton Metre^2) / Kilogram^2 */
  353. #define AU          149504094917.0    /* Metres / Astronomical unit */
  354. #define M_SUN       1.989e30          /* Kilograms / Mass of Sun */
  355. #define YEAR        (365.0 * 24 * 60 * 60) /* Seconds / Year */
  356.  
  357. /*  From Newton's second law, F = ma,
  358.  
  359.            Newton = kg m / sec^2
  360.  
  361.     the fundamental units of the gravitational constant are:
  362.  
  363.            G = N m^2 / kg^2
  364.              = (kg m / sec^2) m^2 / kg^2
  365.              = kg m^3 / sec^2 kg^2
  366.              = m^3 / sec^2 kg
  367.  
  368.     The conversion factor, therefore,  between  the  SI  gravitational
  369.     constant and its equivalent in our units is:
  370.  
  371.            K = AU^3 / YEAR^2 M_SUN
  372.  
  373. */
  374.  
  375. #define GRAV_CONV   ((AU * AU * AU) / ((YEAR * YEAR) * M_SUN))
  376.  
  377. /*  And finally the  gravitational  constant  itself  is  obtained  by
  378.     dividing the SI definition by this conversion factor.  */
  379.  
  380. #define GRAVCON     (G_SI / GRAV_CONV)
  381.  
  382. /*  We also want to come up with approximate sizes for the  objects we
  383.     create.   The  actual  sizes  based  on the density of the objects
  384.     result in everything looking like geometrical points  so,  in  the
  385.     rich  tradition  of  celestial  maps, we enormously exaggerate the
  386.     sizes of objects to render them visible.  Our magic number is  one
  387.     tenth of the cube root of the mass of the object.  */
  388.  
  389. #define DENSCON     0.1
  390.  
  391. static Boolean functional;            /* C:command is returning result */
  392. static ads_real numsteps = 50;        /* Default number of steps to run */
  393. static double simtime = 0.0;          /* Simulated time */
  394. static long stepno = 0;               /* Step number */
  395. static char fhandle[HANDLEN];         /* Handle of reference frame entity */
  396. static int framei = -1;               /* Reference frame object index */
  397.  
  398. /*  Command definition and dispatch table.  */
  399.  
  400. struct {
  401.         char *cmdname;
  402.         void (*cmdfunc)();
  403. } cmdtab[] = {
  404. /*        Name         Function  */
  405. {/*MSG3*/"DEMO",       demo},         /* Create demo case */
  406. {/*MSG4*/"FRAME",      frame},        /* Set local reference frame */
  407. {/*MSG5*/"MASS",       mass},         /* Create new mass */
  408. {/*MSG6*/"MASSEDIT",   massedit},     /* Edit existing mass */
  409. {/*MSG7*/"RESET",      reset},        /* Erase orbital paths */
  410. {/*MSG8*/"RUN",        run},          /* Run simulation */
  411. {/*MSG9*/"SETGRAV",    setgrav},      /* Set mode variables */
  412. {/*MSG10*/"UPDATE",     update}       /* Update masses to last state */
  413. };
  414.  
  415. /*  Particle structure.  */
  416.  
  417. typedef struct {
  418.         char partname[32];            /* Particle name */
  419.         ads_point position;           /* Location in space */
  420.         ads_point velocity;           /* Velocity vector */
  421.         ads_point acceleration;       /* Acceleration vector */
  422.         ads_point lastpos;            /* Last plotted position */
  423.         ads_real mass;                /* Mass */
  424.         ads_real radius;              /* Radius */
  425.         int colour;                   /* Entity colour */
  426.         char partblock[HANDLEN];      /* Block defining particle */
  427. } particle;
  428.  
  429. static particle pt;                   /* Static particle structure */
  430. static particle *ptab = NULL;         /* Particle table */
  431. static int nptab = 0;                 /* Number of in-memory particles */
  432.  
  433. /*  Particle definition block attributes.  */
  434.  
  435. #define ParticleBlock  /*MSG11*/"PARTICLE"  /* Particle block name */
  436. struct {                              /* Attribute tag table */
  437.         char *attname;                /* Attribute tag name */
  438.         ads_real *attvar;             /* Variable address */
  439.         char *attprompt;              /* Prompt for attribute */
  440. } partatt[] = {
  441.         {/*MSG12*/"VELOCITY_X", &pt.velocity[X], /*MSG13*/"Velocity (X)"},
  442.         {/*MSG14*/"VELOCITY_Y", &pt.velocity[Y], /*MSG15*/"Velocity (Y)"},
  443.         {/*MSG16*/"VELOCITY_Z", &pt.velocity[Z], /*MSG17*/"Velocity (Z)"},
  444.         {/*MSG18*/"MASS",       &pt.mass,        /*MSG19*/"Mass"}
  445.        };
  446.  
  447. /*  Modal  variables.   Default   initial   values   below   are   for
  448.     documentation  only.   The actual defaults are reset in initacad()
  449.     upon entry to the drawing editor, so that they will be placed with
  450.     the   default  values  in  the  modal  variable  block  if  it  is
  451.     subsequently created for a new drawing.  */
  452.  
  453. static Boolean drawonly = True,       /* DRAWONLY: Draw, but don't add entities
  454.                                                    if 1.  Make LINEs if 0. */
  455.                showstep = True,       /* SHOWSTEP: Show step in status line. */
  456.                showtime = True;       /* SHOWTIME: Show time in status line. */
  457. static ads_real stepsize = 0.1,       /* STEPSIZE: Motion step size factor. */
  458.                 stepmin = 0.00001;    /* STEPMIN:  Smallest step to use. */
  459.  
  460. /*  Modal attribute definition.  */
  461.  
  462. #define ModalBlock  /*MSG20*/"GRAVITY_MODES" /* Mode variable block name */
  463. struct {
  464.         char *attname;                /* Attribute tag name */
  465.         int attvari;                  /* Variable index */
  466.         char *attprompt;              /* Prompt for variable */
  467. } varatt[] = {
  468.         {/*MSG21*/"DRAWONLY", 1, /*MSG22*/"Output to display only? "},
  469.         {/*MSG23*/"SHOWSTEP", 3, /*MSG24*/"Display step number?    "},
  470.         {/*MSG25*/"SHOWTIME", 2, /*MSG26*/"Display time?           "},
  471.         {/*MSG27*/"STEPSIZE", 4, /*MSG28*/"Step size?              "},
  472.         {/*MSG29*/"STEPMIN",  5, /*MSG30*/"Minimum step (years)?   "}
  473.        };
  474.  
  475. #if 0
  476. /*-----------------------------------------------------------------------*/
  477. /* MAIN -- the main routine */
  478.  
  479. void main(argc, argv)
  480.   int argc;
  481.   char *argv[];
  482. {
  483.     int stat, cindex, scode = RSRSLT;
  484.     char errmsg[80];
  485.  
  486.     ads_init(argc, argv);             /* Initialise the application */
  487.  
  488.     /* Main dispatch loop. */
  489.  
  490.     while (True) {
  491.  
  492.         if ((stat = ads_link(scode)) < 0) {
  493.             sprintf(errmsg,
  494.                     /*MSG31*/"GRAVITY: bad status from ads_link() = %d\n",
  495.                     stat);
  496. #ifdef Macintosh
  497.             macalert(errmsg);
  498. #else
  499.             puts(errmsg);
  500.             fflush(stdout);
  501. #endif /* Macintosh */
  502.             exit(1);
  503.         }
  504.  
  505.         scode = RSRSLT;               /* Default return code */
  506.  
  507.         switch (stat) {
  508.  
  509.         case RQXLOAD:                 /* Load functions.  Called at the start
  510.                                          of the drawing editor.  Re-initialise
  511.                                          the application here. */
  512.             scode = -(funcload() ? RSRSLT : RSERR);
  513.             break;
  514.  
  515.         case RQXUNLD:                 /* Application unload request. */
  516.             break;
  517.  
  518.         case RQSUBR:                  /* Evaluate external lisp function */
  519.             cindex = ads_getfuncode();
  520.             functional = False;
  521.             if (!initacad(False)) {
  522.                 ads_printf(/*MSG32*/"\nUnable to initialise application.\n");
  523.             } else {
  524.  
  525.                 /* Execute the command from the command table with
  526.                    the index associated with this function. */
  527.  
  528.                 if (cindex > 0) {
  529.                     cindex--;
  530.                     assert(cindex < ELEMENTS(cmdtab));
  531.                     (*cmdtab[cindex].cmdfunc)();
  532.                 }
  533.             }
  534.             if (!functional)
  535.                 ads_retvoid();        /* Void result */
  536.             break;
  537.  
  538.         default:
  539.             break;
  540.         }
  541.     }
  542. }
  543. #endif
  544.  
  545. AcRx::AppRetCode acrxEntryPoint(AcRx::AppMsgCode msg, void* ptr)
  546. {
  547.     int stat, cindex, scode = RSRSLT;
  548.     char errmsg[80];
  549.  
  550.  
  551.     if (ptr != NULL) {
  552.         // We have been handed some kind of object
  553.         // but we aren't going to do anything with it.
  554.     }
  555.  
  556.     switch(msg) {
  557.     case AcRx::kInitAppMsg:
  558.         break;
  559.         case AcRx::kInvkSubrMsg:
  560.             cindex = ads_getfuncode();
  561.             functional = False;
  562.             if (!initacad(False)) {
  563.                 ads_printf(/*MSG32*/"\nUnable to initialise application.\n");
  564.             } else {
  565.  
  566.                 /* Execute the command from the command table with
  567.                    the index associated with this function. */
  568.  
  569.                 if (cindex > 0) {
  570.                     cindex--;
  571.                     assert(cindex < ELEMENTS(cmdtab));
  572.                     (*cmdtab[cindex].cmdfunc)();
  573.                 }
  574.             }
  575.             if (!functional)
  576.                 ads_retvoid();        /* Void result */
  577.             break;
  578.         case AcRx::kLoadADSMsg:
  579.                       /* Load functions.  Called at the start
  580.                                          of the drawing editor.  Re-initialise
  581.                                          the application here. */
  582.             scode = -(funcload() ? RSRSLT : RSERR);
  583.             break;
  584.         case AcRx::kUnloadADSMsg:
  585.             ads_printf(/*MSG2*/"Unloading.\n");
  586.             break;
  587.     case AcRx::kUnloadAppMsg:
  588.         default:
  589.         break;
  590.     }
  591.     return AcRx::kRetOK;
  592. }
  593.  
  594.  
  595. /* FUNCLOAD  --  Load external functions into AutoLISP */
  596.  
  597. static Boolean funcload()
  598. {
  599.     char ccbuf[40];
  600.     int i;
  601.  
  602.     strcpy(ccbuf, /*MSG0*/"C:");
  603.     for (i = 0; i < ELEMENTS(cmdtab); i++) {
  604.         strcpy(ccbuf + 2, cmdtab[i].cmdname);
  605.         ads_defun(ccbuf, i + 1);
  606.     }
  607.  
  608.     return initacad(True);            /* Reset AutoCAD initialisation */
  609. }
  610.  
  611. /*  ALLOC  --  Allocate storage and fail if it can't be obtained.  */
  612.  
  613. static char *alloc(unsigned nbytes)
  614. {
  615.     char *cp;
  616.  
  617.     if ((cp = (char *)malloc(nbytes)) == NULL) {
  618. #if 0
  619.         ads_abort(/*MSG33*/"Out of memory");
  620. #endif
  621.  
  622. ;
  623.     }
  624.     return cp;
  625. }
  626.  
  627. /*  DEFMASSBLK  --  Create the mass definition block. */
  628.  
  629. static void defmassblk()
  630. {
  631.     ads_point centre, ucentre;
  632.     ads_real radius = 1;
  633.     ads_point ax, ax1;
  634.     struct resbuf rb1, rb2, ocolour;
  635.     ads_name e1, e2;
  636.     int i;
  637.     ads_name matbss;
  638.     ads_name ename;
  639.     ads_point atx;
  640.  
  641.     Spoint(ucentre, 4, 4, 0);
  642.     rb1.restype = rb2.restype = RTSHORT;
  643.     rb1.resval.rint = 1;              /* From UCS */
  644.     rb2.resval.rint = 0;              /* To world */
  645.     ads_trans(ucentre, &rb1, &rb2, False, centre);
  646.     CommandB();
  647.     ads_command(RTSTR, /*MSG0*/"_.UCS",
  648.                 RTSTR, /*MSG0*/"_X", RTREAL, 90.0, RTNONE);
  649.  
  650.     ads_getvar(/*MSG0*/"CECOLOR", &ocolour);
  651.     /* AutoCAD  currently returns  "human readable" colour strings
  652.        like "1 (red)" for the standard colours.  Trim  the  string
  653.        at  the  first space to guarantee we have a valid string to
  654.        restore the colour later.  */
  655.     if (strchr(ocolour.resval.rstring, ' ') != NULL)
  656.         *strchr(ocolour.resval.rstring, ' ') = EOS;
  657.  
  658.     ads_command(RTSTR, /*MSG0*/"_.COLOUR", RTSTR, /*MSG96*/"_BYBLOCK", RTNONE);
  659.     rb1.resval.rint = 0;              /* From world */
  660.     rb2.resval.rint = 1;              /* To new UCS */
  661.     ads_trans(centre, &rb1, &rb2, False, centre);
  662.     ax[X] = ax1[X] = centre[X];
  663.     ax[Y] = centre[Y] + radius;
  664.     ax[Z] = ax1[Z] = centre[Z];
  665.     ax1[Y] = centre[Y] - radius;
  666.     ads_command(RTSTR, /*MSG0*/"_.LINE", RT3DPOINT, ax, RT3DPOINT, ax1,
  667.                 RTSTR, "", RTNONE);
  668.     ads_entlast(e1);
  669.     ads_command(RTSTR, /*MSG0*/"_.Arc", RT3DPOINT, ax, RTSTR, /*MSG0*/"_E",
  670.                 RT3DPOINT, ax1, RTSTR, /*MSG0*/"_A", RTSTR, "<<180.0",
  671.                 RTNONE);
  672.     ads_entlast(e2);
  673.     PushVar(/*MSG0*/"SURFTAB1", stab1, SphereLong, rint);
  674.     PushVar(/*MSG0*/"SURFTAB2", stab2, SphereLat, rint);
  675.     ads_command(RTSTR, /*MSG0*/"_.REVSURF", RTLB, RTENAME, e2, RT3DPOINT, ax,
  676.                 RTLE, RTLB, RTENAME, e1, RT3DPOINT, centre, RTLE, RTSTR, "",
  677.                 RTSTR, "", RTNONE);
  678.     PopVar(/*MSG0*/"SURFTAB2", stab2);
  679.     PopVar(/*MSG0*/"SURFTAB1", stab1);
  680.     ads_command(RTSTR, /*MSG0*/"_.COLOUR",
  681.                 RTSTR, ocolour.resval.rstring, RTNONE);
  682.     free(ocolour.resval.rstring);
  683.     ads_entdel(e1);
  684.     ads_entdel(e2);
  685.     ads_entlast(e2);
  686.  
  687.     /* Create attributes */
  688.  
  689.     Spoint(atx, 4, 2.75, 0);
  690.     ads_ssadd(e2, NULL, matbss);
  691.  
  692.     PushVar(/*MSG0*/"AFLAGS", saflags, 1, rint);  /* Invisible */
  693.     ads_command(RTSTR, /*MSG0*/"_.ATTDEF", RTSTR, "",
  694.                 RTSTR, /*MSG34*/"PARTNAME",
  695.                 RTSTR, /*MSG35*/"Particle name", RTSTR, "", RT3DPOINT, atx,
  696.                 RTREAL, 0.2, RTREAL, 0.0, RTNONE);
  697.     ads_entlast(ename);
  698.     ads_ssadd(ename, matbss, matbss);
  699.  
  700.     for (i = 0; i < ELEMENTS(partatt); i++) {
  701.         atx[Y] -= 0.25;
  702.         ads_command(RTSTR, /*MSG0*/"_.ATTDEF",
  703.                     RTSTR, "", RTSTR, partatt[i].attname,
  704.                     RTSTR, partatt[i].attprompt, RTSTR, "0", RT3DPOINT, atx,
  705.                     RTREAL, 0.2, RTREAL, 0.0, RTNONE);
  706.         ads_entlast(ename);
  707.         ads_ssadd(ename, matbss, matbss);
  708.     }
  709.     PopVar(/*MSG0*/"AFLAGS", saflags);
  710.  
  711.     /* Collect sphere and attributes into a block. */
  712.  
  713.     ads_command(RTSTR, /*MSG0*/"_.BLOCK",
  714.                 RTSTR, ParticleBlock, RT3DPOINT, centre,
  715.                 RTPICKS, matbss, RTSTR, "", RTNONE);
  716.     ads_command(RTSTR, /*MSG0*/"_.UCS", RTSTR, /*MSG0*/"_Prev", RTNONE);
  717.     CommandE();
  718.     ads_ssfree(matbss);
  719. }
  720.  
  721. /*  ENTITEM  --  Search an entity buffer chain and return an item
  722.                  with the specified group code.  */
  723.  
  724. static struct resbuf *entitem(struct resbuf *rchain, int gcode)
  725. {
  726.     while ((rchain != NULL) && (rchain->restype != gcode))
  727.         rchain = rchain->rbnext;
  728.  
  729.     return rchain;
  730. }
  731.  
  732. /*  ENTINFO  --  Obtain information about a particle entity:
  733.  
  734.                        Handle
  735.                        Position
  736.                        Size (from scale factor of unit block)
  737.                        Colour
  738. */
  739.  
  740. static void entinfo(ads_name ename, char *h, ads_point p, ads_real *r, int *c)
  741. {
  742.     struct resbuf *rent, *rh;
  743.  
  744.     rent = ads_entget(ename);
  745.     if ((rh = entitem(rent, 5)) != NULL)
  746.         strcpy(h, rh->resval.rstring);
  747.     else
  748.         h[0] = EOS;
  749.     rh = entitem(rent, 10);
  750.     assert(rh != NULL);
  751.     Cpoint(p, rh->resval.rpoint);
  752.     rh = entitem(rent, 41);
  753.     assert(rh != NULL);
  754.     *r = rh->resval.rreal;
  755.     if ((rh = entitem(rent, 62)) != NULL) {
  756.         *c = rh->resval.rint;
  757.         if (*c == 0)                  /* Naked colour by block? */
  758.             *c = 7;                   /* Forbidden: make it white.  Q.C.D. */
  759.     } else {
  760.         /* This entity derives its colour from the layer.  Get
  761.            the colour from the layer table. */
  762.         *c = 7;
  763.         if ((rh = entitem(rent, 8)) != NULL) {
  764.             if ((rh = ads_tblsearch(/*MSG0*/"LAYER", rh->resval.rstring,
  765.                                     False)) != NULL) {
  766.                 struct resbuf *lh = entitem(rh, 62);
  767.                 if (lh != NULL) {
  768.                     int lc = abs(lh->resval.rint);
  769.                     if (lc > 0 && lc < 256) {
  770.                         *c = lc;
  771.                     }
  772.                 }
  773.                 ads_relrb(rh);
  774.             }
  775.         }
  776.     }
  777.     ads_relrb(rent);
  778. }
  779.  
  780. /*  PARTEXT  --  Extract particles present in drawing and build the
  781.                  in-memory particle table.  */
  782.  
  783. static void partext()
  784. {
  785.     long i, l;
  786.     struct resbuf rbet, rbb;
  787.     struct resbuf *rb, *rent, *vb;
  788.     ads_name ename, vname;
  789.  
  790.     if (ptab != NULL) {
  791.         free(ptab);
  792.     }
  793.     ptab = NULL;
  794.     nptab = 0;
  795.  
  796.     /* Build the SSGET entity buffer chain to filter for block
  797.        insertions of the named block on the named layer. */
  798.  
  799.     rbet.restype = 0;                 /* Entity type */
  800.     rbet.resval.rstring = /*MSG0*/"INSERT";
  801.     rbet.rbnext = &rbb;
  802.     rbb.restype = 2;                  /* Block name */
  803.     rbb.resval.rstring = ParticleBlock;
  804.     rbb.rbnext = NULL;
  805.  
  806.     if (ads_ssget(/*MSG0*/"_X", NULL, NULL, &rbet, vname) != RTNORM) {
  807.         return;                       /* No definitions in database */
  808.     }
  809.  
  810.     /* We found one or more definitions.  Process  the  attributes
  811.        that  follow  each, plugging the values into material items
  812.        which get attached to the in-memory definition chain.  */
  813.  
  814.     if (ads_sslength(vname, &l) < 0)
  815.         l = 0;
  816.  
  817.     nptab = l;                        /* Save particle count */
  818.     ptab = (particle *) alloc(nptab * sizeof(particle));
  819.     for (i = 0; i < l; i++) {
  820.         ads_name pname;
  821.  
  822.         ads_ssname(vname, i, ename);
  823.         memcpy(pname, ename, sizeof ename);
  824.  
  825.         memset(&pt, 0, sizeof pt);
  826.         while (True) {
  827.             ads_entnext(ename, ename);
  828.             rent = ads_entget(ename);
  829.             if (rent == NULL) {
  830.                 ads_printf(/*MSG36*/"PARTEXT: Can't read attribute.\n");
  831.                 break;
  832.             }
  833.             rb = entitem(rent, 0);    /* Entity type */
  834.             if (rb != NULL) {
  835.                 if (strcmp(rb->resval.rstring, /*MSG0*/"ATTRIB") != 0)
  836.                     break;
  837.                 rb = entitem(rent, 2);  /* Attribute tag */
  838.                 vb = entitem(rent, 1);  /* Attribute value */
  839.                 if (rb != NULL && vb != NULL) {
  840.                     if (strcmp(rb->resval.rstring, /*MSG37*/"PARTNAME") == 0) {
  841.                         strcpy(pt.partname, vb->resval.rstring);
  842.                     } else {
  843.                         int j;
  844.  
  845.                         for (j = 0; j < ELEMENTS(partatt); j++) {
  846.                             if (strcmp(rb->resval.rstring,
  847.                                        partatt[j].attname) == 0) {
  848.                                 *partatt[j].attvar = atof(vb->resval.rstring);
  849.                                 break;
  850.                             }
  851.                         }
  852.                     }
  853.                 }
  854.             }
  855.             ads_relrb(rent);
  856.         }
  857.         entinfo(pname, pt.partblock, pt.position, &pt.radius, &pt.colour);
  858.         memcpy(&(ptab[(int) i]), &pt, sizeof(particle));
  859.     }
  860.     ads_ssfree(vname);                /* Release selection set */
  861. }
  862.  
  863. /*  ADDVEC  --  Add two vectors, a = b + c  */
  864.  
  865. static void addvec(ads_real *ap, ads_real *bp, ads_real *cp)
  866. {
  867.     ap[X] = bp[X] + cp[X];
  868.     ap[Y] = bp[Y] + cp[Y];
  869.     ap[Z] = bp[Z] + cp[Z];
  870. }
  871.  
  872. /*  SUBVEC  --  Subtract two vectors, a = b - c  */
  873.  
  874. static void subvec(ads_real *ap, ads_real *bp, ads_real *cp)
  875. {
  876.     ap[X] = bp[X] - cp[X];
  877.     ap[Y] = bp[Y] - cp[Y];
  878.     ap[Z] = bp[Z] - cp[Z];
  879. }
  880.  
  881. /*  SQABSV  --  Square of absolute value of a vector.  */
  882.  
  883. static ads_real sqabsv(ads_real *ap)
  884. {
  885.     return (ap[X] * ap[X] + ap[Y] * ap[Y] + ap[Z] * ap[Z]);
  886. }
  887.  
  888. /*  SUMVEC  --  Add a linear multiple to another vector, a = b + t * c.  */
  889.  
  890. static void sumvec(ads_real *ap, ads_real *bp, ads_real t, ads_real *cp)
  891. {
  892.     ap[X] = bp[X] + t * cp[X];
  893.     ap[Y] = bp[Y] + t * cp[Y];
  894.     ap[Z] = bp[Z] + t * cp[Z];
  895. }
  896.  
  897. /*  VARBLOCKDEF  --  Create the block that carries the attributes that
  898.                      define the modal variables. */
  899.  
  900. static void varblockdef()
  901. {
  902.     int i;
  903.     ads_name varbss;
  904.     ads_name ename;
  905.     ads_point atx;
  906.  
  907.     atx[X] = 4.0;
  908.     atx[Y] = 4.0;
  909.     atx[Z] = 0.0;
  910.     ads_ssadd(NULL, NULL, varbss);
  911.  
  912.     CommandB();
  913.     PushVar(/*MSG0*/"AFLAGS", saflags, 1, rint);  /* Invisible */
  914.  
  915.     for (i = 0; i < ELEMENTS(varatt); i++) {
  916.         atx[Y] -= 0.25;
  917.         ads_command(RTSTR, /*MSG0*/"_.ATTDEF",
  918.                     RTSTR, "", RTSTR, varatt[i].attname,
  919.                     RTSTR, varatt[i].attprompt, RTSTR, "0", RT3DPOINT, atx,
  920.                     RTREAL, 0.2, RTREAL, 0.0, RTNONE);
  921.         ads_entlast(ename);
  922.         ads_ssadd(ename, varbss, varbss);
  923.     }
  924.     PopVar(/*MSG0*/"AFLAGS", saflags);
  925.  
  926.     ads_command(RTSTR, /*MSG0*/"_.BLOCK", RTSTR, ModalBlock, RTSTR, "4,4",
  927.                 RTPICKS, varbss, RTSTR, "", RTNONE);
  928.     CommandE();
  929.     ads_ssfree(varbss);
  930. }
  931.  
  932. /*  SAVEMODES  --  Save application modes as attributes of the mode
  933.                    block.  */
  934.  
  935. static void savemodes()
  936. {
  937.     int i;
  938.     struct resbuf rbb;
  939.     ads_name ename, vname;
  940.     long l;
  941.  
  942.     /* Build the SSGET entity buffer chain to filter for block
  943.        insertions of the named block on the named layer. */
  944.  
  945.     rbb.restype = 2;                  /* Block name */
  946.     rbb.resval.rstring = ModalBlock;
  947.     rbb.rbnext = NULL;
  948.  
  949.     if (ads_ssget(/*MSG0*/"_X", NULL, NULL, &rbb, vname) == RTNORM) {
  950.  
  951.         /* Delete all definitions found. */
  952.  
  953.         if (ads_sslength(vname, &l) < 0)
  954.             l = 0;
  955.  
  956.         if (l > 0) {
  957.             CommandB();
  958.             ads_command(RTSTR, /*MSG0*/"_.ERASE",
  959.                         RTPICKS, vname, RTSTR, "", RTNONE);
  960.             CommandE();
  961.         }
  962.         ads_ssfree(vname);            /* Release selection set */
  963.     }
  964.  
  965.     /* Now insert the modal variable block, attaching the
  966.        mode variables to its attributes. */
  967.  
  968.     CommandB();
  969.     ads_command(RTSTR, /*MSG0*/"_.INSERT", RTSTR, ModalBlock,
  970.                 RTSTR, "0,0", RTSTR, "1", RTSTR, "1", RTSTR, "0",
  971.                 RTNONE);
  972.     for (i = 0; i < ELEMENTS(varatt); i++) {
  973.         char attval[20];
  974.  
  975. #define YorN(x) ((x) ? /*MSG38*/"Yes" : /*MSG39*/"No")
  976.         switch (varatt[i].attvari) {
  977.         case 1:
  978.             strcpy(attval, YorN(drawonly));
  979.             break;
  980.         case 2:
  981.             strcpy(attval, YorN(showtime));
  982.             break;
  983.         case 3:
  984.             strcpy(attval, YorN(showstep));
  985.             break;
  986.         case 4:
  987.             sprintf(attval, "%.12g", stepsize);
  988.             break;
  989.         case 5:
  990.             sprintf(attval, "%.12g", stepmin);
  991.             break;
  992.         }
  993. #undef YorN
  994.         ads_command(RTSTR, attval, RTNONE);
  995.     }
  996.     ads_entlast(ename);
  997.     ads_command(RTSTR, /*MSG0*/"_.CHPROP", RTENAME, ename, RTSTR, "",
  998.                 RTSTR, /*MSG0*/"_layer", RTSTR, FrozenLayer, RTSTR, "",
  999.                 RTNONE);
  1000.     CommandE();
  1001. }
  1002.  
  1003. /*  BOOLVAR  --  Determine Boolean value from attribute string.  We
  1004.                  recognise numbers (nonzero means True) and the
  1005.                  strings "True", "False", "Yes", and "No",
  1006.                  in either upper or lower case.  */
  1007.  
  1008. static Boolean boolvar(char *varname, char *s)
  1009. {
  1010.     char ch = *s;
  1011.  
  1012.     if (islower(ch))
  1013.         ch = toupper(ch);
  1014.     if (isdigit(ch)) {
  1015.         return (atoi(s) != 0) ? True : False;
  1016.     }
  1017.  
  1018.     switch (ch) {
  1019.     case /*MSG40*/'Y':
  1020.     case /*MSG41*/'T':
  1021.         return True;
  1022.     case /*MSG42*/'N':
  1023.     case /*MSG43*/'F':
  1024.         return False;
  1025.     default:
  1026.         ads_printf(/*MSG44*/"Invalid Boolean value for %s: %s.\n", varname, s);
  1027.         break;
  1028.     }
  1029.     return False;
  1030. }
  1031.  
  1032. /*  VARSET  --  Set modal variables from the block attributes.  */
  1033.  
  1034. static void varset()
  1035. {
  1036.     struct resbuf rbb;
  1037.     ads_name vname;
  1038.     long l;
  1039.  
  1040.     /* Build the SSGET entity buffer chain to filter for block
  1041.        insertions of the named block on the named layer. */
  1042.  
  1043.     rbb.restype = 2;                  /* Block name */
  1044.     rbb.resval.rstring = ModalBlock;
  1045.     rbb.rbnext = NULL;
  1046.  
  1047.     if (ads_ssget(/*MSG0*/"_X", NULL, NULL, &rbb, vname) == RTNORM) {
  1048.  
  1049.         if (ads_sslength(vname, &l) < 0)
  1050.             l = 0;
  1051.  
  1052.         if (l > 0) {
  1053.             ads_name ename;
  1054.  
  1055.             if (ads_ssname(vname, 0L, ename) == RTNORM) {
  1056.                 while (ads_entnext(ename, ename) == RTNORM) {
  1057.                     int i;
  1058.                     struct resbuf *rb = ads_entget(ename),
  1059.                                   *et, *at, *av;
  1060.  
  1061.                     if (rb == NULL)
  1062.                         break;
  1063.                     et = entitem(rb, 0);
  1064.                     assert(et != NULL);
  1065.                     if (strcmp(et->resval.rstring, /*MSG0*/"ATTRIB") == 0) {
  1066.                         at = entitem(rb, 2);  /* Attribute tag */
  1067.                         av = entitem(rb, 1);  /* Attribute value */
  1068.                         if (at == NULL || av == NULL)
  1069.                             break;
  1070.                         for (i = 0; i < ELEMENTS(varatt); i++) {
  1071.                             if (strcmp(at->resval.rstring,
  1072.                                        varatt[i].attname) == 0) {
  1073.                                 switch (varatt[i].attvari) {
  1074.                                 case 1:
  1075.                                     drawonly = boolvar(at->resval.rstring,
  1076.                                                        av->resval.rstring);
  1077.                                     break;
  1078.                                 case 2:
  1079.                                     showtime = boolvar(at->resval.rstring,
  1080.                                                        av->resval.rstring);
  1081.                                     break;
  1082.                                 case 3:
  1083.                                     showstep = boolvar(at->resval.rstring,
  1084.                                                        av->resval.rstring);
  1085.                                     break;
  1086.                                 case 4:
  1087.                                     stepsize = 0.1;    /* Default if error */
  1088.                                     sscanf(av->resval.rstring, "%lf",
  1089.                                            &stepsize);
  1090.                                     break;
  1091.                                 case 5:
  1092.                                     stepmin = 0.00001; /* Default if error */
  1093.                                     sscanf(av->resval.rstring, "%lf",
  1094.                                            &stepmin);
  1095.                                     break;
  1096.                                 }
  1097.                             }
  1098.                         }
  1099.                     } else if (strcmp(et->resval.rstring,
  1100.                                       /*MSG0*/"SEQEND") == 0) {
  1101.                         ads_relrb(rb);
  1102.                         break;
  1103.                     }
  1104.                     ads_relrb(rb);
  1105.                 }
  1106.             }
  1107.         }
  1108.         ads_ssfree(vname);            /* Release selection set */
  1109.     }
  1110. }
  1111.  
  1112. /*  SETFRAME  --  Set reference frame index from handle of
  1113.                   reference frame entity in effect.  */
  1114.  
  1115. static void setframe()
  1116. {
  1117.     int i;
  1118.  
  1119.     framei = -1;
  1120.     for (i = 0; i < nptab; i++) {
  1121.         if (strcmp(ptab[i].partblock, fhandle) == 0) {
  1122.             framei = i;
  1123.             break;
  1124.         }
  1125.     }
  1126. }
  1127.  
  1128. /*  INITACAD  --  Initialise the required modes in the AutoCAD
  1129.                   drawing.  */
  1130.  
  1131. static Boolean initacad(Boolean reset)
  1132. {
  1133.     static Boolean initdone, initok;
  1134.  
  1135.     if (reset) {
  1136.         initdone = False;
  1137.         initok = True;
  1138.     } else {
  1139.         if (!initdone) {
  1140.             struct resbuf rb;
  1141.             struct resbuf *ep;
  1142.  
  1143.             initok = False;
  1144.  
  1145.             /* Reset the program modes to standard values upon
  1146.                entry to the drawing editor. */
  1147.  
  1148.             stepno = 0;               /* Step 0 */
  1149.             numsteps = 50;            /* Default number of steps */
  1150.             simtime = 0.0;            /* No elapsed time */
  1151.             stepsize = 0.1;           /* Default step size */
  1152.             stepmin = 0.00001;        /* Default minimum time: none */
  1153.             drawonly = True;          /* Draw using ads_grdraw() */
  1154.             showstep = True;          /* Show step number */
  1155.             showtime = True;          /* Show elapsed time */
  1156.             framei = -1;              /* No reference frame object */
  1157.             fhandle[0] = EOS;         /* Clear reference frame handle */
  1158.  
  1159.             /* Enable  handles  if  they aren't  already on.  We use
  1160.                handles to link  from  our  in-memory  table  to  the
  1161.                masses  in  the  AutoCAD  database,  and we also keep
  1162.                track of the reference frame entity  by  its  handle. */
  1163.  
  1164.             ads_getvar(/*MSG0*/"HANDLES", &rb);
  1165.             if (!rb.resval.rint) {
  1166.                 CommandB();
  1167.                 ads_command(RTSTR, /*MSG0*/"_.handles",
  1168.                             RTSTR, /*MSG0*/"_on", RTNONE);
  1169.                 CommandE();
  1170.                 ads_getvar(/*MSG0*/"HANDLES", &rb);
  1171.                 if (!rb.resval.rint) {
  1172.                     ads_printf(/*MSG45*/"Cannot enable handles.\n");
  1173.                     initdone = True;
  1174.                     return (initok = False);
  1175.                 }
  1176.             }
  1177.  
  1178.             /* Create required drawing objects. */
  1179.  
  1180.             /* If there isn't already a "frozen solid" layer, create one. */
  1181.  
  1182.             if ((ep = ads_tblsearch(/*MSG0*/"LAYER",
  1183.                                     FrozenLayer, False)) == NULL) {
  1184.                 CommandB();
  1185.                 ads_command(RTSTR, /*MSG0*/"_.LAYER",
  1186.                             RTSTR, /*MSG0*/"_NEW", RTSTR, FrozenLayer,
  1187.                             RTSTR, /*MSG0*/"_FREEZE",
  1188.                             RTSTR, FrozenLayer, RTSTR, "", RTNONE);
  1189.                 ads_command(RTSTR, /*MSG0*/"_.LAYER",
  1190.                             RTSTR, /*MSG0*/"_NEW", RTSTR, OrbitLayer,
  1191.                             RTSTR, "", RTNONE);
  1192.                 CommandE();
  1193.                 defmassblk();         /* Define the particle block */
  1194.             } else {
  1195.                 ads_relrb(ep);
  1196.             }
  1197.  
  1198.             /* Create the block definition for our modal variables. */
  1199.  
  1200.             if ((ep = ads_tblsearch(/*MSG0*/"BLOCK",
  1201.                                     ModalBlock, False)) == NULL) {
  1202.                 varblockdef();
  1203.                 savemodes();
  1204.             } else {
  1205.                 ads_relrb(ep);
  1206.                 varset();             /* Load modals from block */
  1207.             }
  1208.  
  1209.             initdone = initok = True;
  1210.         }
  1211.     }
  1212.     return initok;
  1213. }
  1214.  
  1215. /*  ADDMASS  --  Insert mass block with attributes.  */
  1216.  
  1217. static void addmass(char *mname, ads_point pos, ads_point vel, ads_real mass)
  1218. {
  1219.     ads_real size;
  1220.     char velx[32], vely[32], velz[32], smas[32];
  1221.  
  1222.     size = cuberoot(mass) * DENSCON;  /* Set size by cube root of mass */
  1223.     CommandB();
  1224.     sprintf(velx, "%.12g", vel[X]);
  1225.     sprintf(vely, "%.12g", vel[Y]);
  1226.     sprintf(velz, "%.12g", vel[Z]);
  1227.     sprintf(smas, "%.12g", mass);
  1228.     ads_command(RTSTR, /*MSG0*/"_.INSERT", RTSTR, ParticleBlock, RTPOINT, pos,
  1229.                 RTREAL, size, RTREAL, size, RTREAL, 0.0,
  1230.                 RTSTR, mname,
  1231.                 RTSTR, velx, RTSTR, vely, RTSTR, velz,
  1232.                 RTSTR, smas, RTNONE);
  1233.     CommandE();
  1234. }
  1235.  
  1236. /*  MASS  --  Add mass to database.  */
  1237.  
  1238. static void mass()
  1239. {
  1240.     int k;
  1241.     ads_point pos, vel;
  1242.     ads_real mass;
  1243.     char pname[134];
  1244.  
  1245.     if (ads_getstring(True, /*MSG46*/"\nParticle name (CR at end): ",
  1246.                       pname) != RTNORM)
  1247.         return;
  1248.     ads_initget(1 + 8 + 16, NULL);
  1249.     if (ads_getpoint(NULL, /*MSG47*/"\nPosition: ", pos) != RTNORM)
  1250.         return;
  1251.     ads_initget(1 + 8 + 16, /*MSG48*/"Vector");
  1252.     switch (ads_getpoint(NULL, /*MSG49*/"\nVelocity (or Vector) in AU/Year: ",
  1253.                          vel)) {
  1254.     case RTNORM:
  1255.         break;
  1256.     case RTKWORD:                     /* Vector keyword */
  1257.         if (ads_getpoint(pos, /*MSG50*/"\nVelocity vector: ", vel) != RTNORM)
  1258.             return;
  1259.         for (k = X; k <= Z; k++)
  1260.             vel[k] -= pos[k];
  1261.         break;
  1262.     default:
  1263.         return;
  1264.     }
  1265.     if (ads_getreal(/*MSG51*/"\nMass in suns: ", &mass) != RTNORM)
  1266.         return;
  1267.  
  1268.     addmass(pname, pos, vel, mass);
  1269. }
  1270.  
  1271. /*  DEMO  --  Load a canned demo world.  */
  1272.  
  1273. static void demo()
  1274. {
  1275.     int i, j;
  1276.  
  1277.     typedef struct {                  /* Demo mass table item */
  1278.        char *dmname;                  /* Mass name */
  1279.        ads_point dmpos;               /* Position */
  1280.        ads_point dmvel;               /* Velocity */
  1281.        ads_real dmmass;               /* Mass */
  1282.        int dmcolour;                  /* Colour */
  1283.     } dmass;
  1284.  
  1285.     typedef struct {
  1286.        char *dname;                   /* Name of demo case */
  1287.        ads_point dwind1,              /* Display window corners */
  1288.                  dwind2;
  1289.        ads_real dnstep;               /* Default number of steps */
  1290.        ads_real dssize;               /* Step size */
  1291.        ads_real dssmin;               /* Minimum step size */
  1292.        char *dframe;                  /* Reference frame */
  1293.        dmass *dmtab;                  /* Table of masses */
  1294.        int dmtabl;                    /* Number of masses in table */
  1295.     } demodesc;
  1296.  
  1297.     static dmass interlom[] = {
  1298.        {/*MSG52*/"Interloper", {-167.5, 168, 0}, {2,-1.5,0}, 2, 2},
  1299.        {/*MSG53*/"Star 3", {100, 50, 0}, {-1, 0, 0}, 20, 3},
  1300.        {/*MSG54*/"Star 2", {0, 100, 0}, {2, 0, 0}, 10, 1},
  1301.        {/*MSG55*/"Star 1", {0, 0, 0}, {-3, 0, 0}, 10, 5}
  1302.       };
  1303.     static dmass solarsm[] = {
  1304.        {/*MSG56*/"Sun", {0, 0, 0}, {0, 0, 0}, 1, 2},
  1305.        {/*MSG57*/"Mercury", {0.387, 0, 0}, {0, 10.0965, 0}, 1.666667e-7, 7},
  1306.        {/*MSG58*/"Venus", {0.723, 0, 0}, {0, 7.38628, 0}, 2.44786e-6, 7},
  1307.        {/*MSG59*/"Earth", {1, 0, 0}, {0, 6.21318531, 0}, 3.04044e-6, 5},
  1308.        {/*MSG60*/"Moon", {1,-0.00256952,0}, {0.215831,6.21318531,0},
  1309.         3.6944e-8, 7},
  1310.        {/*MSG61*/"Mars", {1.524, 0, 0}, {0, 5.0894, 0}, 3.22716e-7, 1},
  1311.        {/*MSG62*/"Jupiter", {5.203, 0, 0}, {0, 2.75456, 0}, 0.000954782, 7},
  1312.        {/*MSG63*/"Saturn", {9.539, 0, 0}, {0, 2.03534, 0}, 0.00285837, 2},
  1313.        {/*MSG64*/"Uranus", {19.182, 0, 0}, {0, 1.43423, 0}, 4.38596e-5, 4},
  1314.        {/*MSG65*/"Neptune", {30.058, 0, 0}, {0, 1.14527, 0}, 5.18135e-5, 4}
  1315.       };
  1316.     static dmass nemesis = {
  1317.        /*MSG66*/"Nemesis", {28.8879, 1.67301, 0}, {-65, 0, 0}, 8, 2
  1318.       };
  1319.     static demodesc dmt[] = {
  1320.        {/*MSG67*/"Interloper", {-176.3, -74, 0}, {176, 180, 0},
  1321.         3000, 0.1, 0.00001, NULL,
  1322.         interlom, ELEMENTS(interlom)},
  1323.        {/*MSG68*/"Solar system", {-3.35, -11.82, 0}, {31.9, 13.6, 0},
  1324.         -1, 100, 0, NULL,
  1325.         solarsm, ELEMENTS(solarsm)},
  1326.        {/*MSG69*/"Nemesis", {-1.59, -1.3, 0}, {2.67, 1.79, 0},
  1327.         -1, 100, 0, /*MSG70*/"Sun",
  1328.         solarsm, ELEMENTS(solarsm)}
  1329.       };
  1330.  
  1331.     if (nptab > 0) {
  1332.         ads_printf(/*MSG71*/"\n\
  1333. Demo can be created only in an empty drawing.\n");
  1334.         return;
  1335.     }
  1336.  
  1337.     partext();
  1338.     ads_textscr();                    /* Flip to text screen */
  1339.     ads_printf(/*MSG72*/"Available demonstration cases:\n");
  1340.     for (i = 0; i < ELEMENTS(dmt); i++) {
  1341.         ads_printf("\n%d.  %s", i + 1, dmt[i].dname);
  1342.     }
  1343.     ads_printf("\n");
  1344.     ads_initget(2 + 4, NULL);
  1345.     switch (ads_getint(/*MSG73*/"\nWhich demonstration? ", &i)) {
  1346.     case RTNORM:
  1347.         i--;
  1348.         if (i < ELEMENTS(dmt))
  1349.             break;
  1350.     default:
  1351.         return;
  1352.     }
  1353.  
  1354.     CommandB();
  1355.     ads_command(RTSTR, /*MSG0*/"_.ZOOM", RTSTR, /*MSG0*/"_Window",
  1356.                 RTPOINT, dmt[i].dwind1, RTPOINT, dmt[i].dwind2, RTNONE);
  1357.     CommandE();
  1358.     for (j = 0; j < dmt[i].dmtabl; j++) {
  1359.         CommandB();
  1360.         ads_command(RTSTR, /*MSG0*/"_.COLOUR",
  1361.                     RTSHORT, dmt[i].dmtab[j].dmcolour, RTNONE);
  1362.         CommandE();
  1363.         addmass(dmt[i].dmtab[j].dmname, dmt[i].dmtab[j].dmpos,
  1364.                 dmt[i].dmtab[j].dmvel, dmt[i].dmtab[j].dmmass);
  1365.     }
  1366.  
  1367.     /* Special case for Nemesis demo.  Load the standard Solar
  1368.        System definitions above, then jam the Nemesis mass
  1369.        on top of it.  Here comes trouble.... */
  1370.     if (i == 2) {
  1371.         CommandB();
  1372.         ads_command(RTSTR, /*MSG0*/"_.COLOUR",
  1373.                     RTSHORT, nemesis.dmcolour, RTNONE);
  1374.         CommandE();
  1375.         addmass(nemesis.dmname, nemesis.dmpos, nemesis.dmvel,
  1376.                 nemesis.dmmass);
  1377.     }
  1378.     /* We now return to our elegant program, already in progress. */
  1379.  
  1380.     numsteps = dmt[i].dnstep;
  1381.     stepsize = dmt[i].dssize;
  1382.     stepmin = dmt[i].dssmin;
  1383.     savemodes();
  1384.     CommandB();
  1385.     ads_command(RTSTR, /*MSG0*/"_.COLOUR", RTSTR, /*MSG127*/"_BYLAYER", RTNONE);
  1386.     CommandE();
  1387.     partext();
  1388.  
  1389.     /* If this demo requests a default reference frame, place it
  1390.        in effect. */
  1391.  
  1392.     if (dmt[i].dframe != NULL) {
  1393.         for (j = 0; j < nptab; j++) {
  1394.             if (strcmp(ptab[j].partname, dmt[i].dframe) == 0) {
  1395.                 strcpy(fhandle, ptab[j].partblock);
  1396.                 framei = j;
  1397.                 break;
  1398.             }
  1399.         }
  1400.         assert(j < nptab);
  1401.     }
  1402. }
  1403.  
  1404. /*  MASSEDIT  --  Edit properties of an existing mass.  */
  1405.  
  1406. static void massedit()
  1407. {
  1408.     int i = nptab + 1;
  1409.     ads_name ename;
  1410.     ads_point ppt;
  1411.  
  1412.     partext();                        /* Update in-memory mass table */
  1413.     if (ads_entsel(/*MSG74*/"Select mass to edit: ", ename, ppt) == RTNORM) {
  1414.         struct resbuf *eb = ads_entget(ename), *ha;
  1415.  
  1416.         if ((ha = entitem(eb, 5)) != NULL) {
  1417.             for (i = 0; i < nptab; i++) {
  1418.                 if (strcmp(ptab[i].partblock, ha->resval.rstring) == 0) {
  1419.                     break;
  1420.                 }
  1421.             }
  1422.         }
  1423.     }
  1424.     if (i < nptab) {
  1425.         if (stepno > 0) {
  1426.             reset();
  1427.         }
  1428.         CommandB();
  1429.         ads_command(RTSTR, /*MSG0*/"_.DDATTE", RTENAME, ename, RTNONE);
  1430.         CommandE();
  1431.         partext();
  1432.     } else {
  1433.         ads_printf(/*MSG75*/"\nNo mass selected.\n");
  1434.     }
  1435. }
  1436.  
  1437. /*  RESET  --  Reset simulation, erasing all orbital paths.  */
  1438.  
  1439. static void reset()
  1440. {
  1441.     struct resbuf rbb;
  1442.     ads_name vname;
  1443.     long l;
  1444.  
  1445.     /* Build the SSGET entity buffer chain to select all
  1446.        entities on the orbital path layer. */
  1447.  
  1448.     rbb.restype = 8;                  /* Layer name */
  1449.     rbb.resval.rstring = OrbitLayer;
  1450.     rbb.rbnext = NULL;
  1451.  
  1452.     if (ads_ssget(/*MSG0*/"_X", NULL, NULL, &rbb, vname) == RTNORM) {
  1453.  
  1454.         /* We found one or more definitions.  Delete them. */
  1455.  
  1456.         if (ads_sslength(vname, &l) < 0)
  1457.             l = 0;
  1458.  
  1459.         if (l > 0) {
  1460.             CommandB();
  1461.             ads_command(RTSTR, /*MSG0*/"_.ERASE",
  1462.                         RTPICKS, vname, RTSTR, "", RTNONE);
  1463.             CommandE();
  1464.         }
  1465.         ads_ssfree(vname);            /* Release selection set */
  1466.     }
  1467.     if (drawonly) {
  1468.         ads_grtext(-3, NULL, 0);      /* Clear time display */
  1469.         ads_redraw(NULL, 0);
  1470.     }
  1471.     stepno = 0;                       /* Reset step number */
  1472. }
  1473.  
  1474. /*  FRAME  --  Set reference frame for motion display.  */
  1475.  
  1476. static void frame()
  1477. {
  1478.     ads_name ename;
  1479.     ads_point ppt;
  1480.  
  1481.     fhandle[0] = EOS;                 /* Assume inertial frame */
  1482.     if (ads_entsel(/*MSG76*/"\
  1483. Select reference frame mass or RETURN for inertial: ",
  1484.                    ename, ppt) == RTNORM) {
  1485.         struct resbuf *eb = ads_entget(ename), *ha;
  1486.  
  1487.         if ((ha = entitem(eb, 5)) != NULL) {
  1488.             strcpy(fhandle, ha->resval.rstring);
  1489.         }
  1490.     }
  1491.     setframe();
  1492.     ads_printf(/*MSG77*/"\nReference frame: %s\n",
  1493.                framei < 0 ? /*MSG78*/"Inertial" :
  1494.                 ptab[framei].partname);
  1495. }
  1496.  
  1497. /*  RUN  --  Run simulation.  This can be called as an interactive
  1498.              command, RUN, or functionally as (C:RUN <steps/-time>).  */
  1499.  
  1500. static void run()
  1501. {
  1502.     int i, j, k, n = 0, lastcolour = -1;
  1503.     Boolean resume = (stepno > 0) ? True : False, timelimit;
  1504.     ads_real r, deltat, endtime = 0;
  1505.     char rpt[80];
  1506.     struct resbuf clayer, ocolour;
  1507.     struct resbuf *rb = ads_getargs();
  1508.  
  1509.     /* If an argument was supplied, set the simulation length
  1510.        from it rather than asking the user interactively. */
  1511.  
  1512.     if (rb != NULL) {
  1513.         Boolean argok = True;
  1514.  
  1515.         switch (rb->restype) {
  1516.         case RTREAL:
  1517.             numsteps = rb->resval.rreal;
  1518.             break;
  1519.         case RTSHORT:
  1520.             numsteps = rb->resval.rint;
  1521.             break;
  1522.         default:
  1523.             ads_fail(/*MSG79*/"RUN: Improper argument type");
  1524.             argok = False;
  1525.             break;
  1526.         }
  1527.         if (argok) {
  1528.             if (numsteps == 0) {
  1529.                 ads_fail(/*MSG80*/"RUN: Argument must be nonzero");
  1530.                 argok = False;
  1531.             } else {
  1532.                 if (rb->rbnext != NULL) {
  1533.                     ads_fail(/*MSG81*/"RUN: Too many arguments");
  1534.                     argok = False;
  1535.                 }
  1536.             }
  1537.         }
  1538.         if (!argok)
  1539.             return;
  1540.         functional = True;            /* Mark called as a function */
  1541.     } else {
  1542.  
  1543.         /* Ask user for length of simulation in years or number of
  1544.            integration steps to perform. */
  1545.  
  1546.         sprintf(rpt, /*MSG82*/"Steps or (-simulated years) <%.12g>: ",
  1547.                 numsteps);
  1548.         ads_initget(2, NULL);
  1549.         switch (ads_getreal(rpt, &r)) {
  1550.         case RTCAN:                   /* Control C */
  1551.             return;
  1552.         case RTNONE:                  /* Null input */
  1553.             break;
  1554.         case RTNORM:                  /* Number entered */
  1555.             numsteps = r;
  1556.             break;
  1557.         }
  1558.     }
  1559.  
  1560.     /* If we're starting a new simulation rather than resuming one
  1561.        already underway, reset  the  simulated  time  counter  and
  1562.        extract the current particle information from the database. */
  1563.  
  1564.     if (!resume) {
  1565.         simtime = 0.0;
  1566.         partext();
  1567.     }
  1568.  
  1569.     setframe();                       /* Activate reference frame */
  1570.     if (numsteps > 0) {
  1571.         timelimit = False;
  1572.         n = numsteps;
  1573.     } else {
  1574.         timelimit = True;
  1575.         endtime = simtime - numsteps;
  1576.     }
  1577.  
  1578.     /* Save original layer and colour */
  1579.  
  1580.     ads_getvar(/*MSG0*/"CLAYER", &clayer);
  1581.     ads_getvar(/*MSG0*/"CECOLOR", &ocolour);
  1582.     /* AutoCAD  currently returns  "human readable" colour strings
  1583.        like "1 (red)" for the standard colours.  Trim  the  string
  1584.        at  the  first space to guarantee we have a valid string to
  1585.        restore the colour later.  */
  1586.     if (strchr(ocolour.resval.rstring, ' ') != NULL)
  1587.         *strchr(ocolour.resval.rstring, ' ') = EOS;
  1588.  
  1589.     CommandB();
  1590.     ads_command(RTSTR, /*MSG0*/"_.UCS", RTSTR, /*MSG0*/"_World", RTNONE);
  1591.     if (!drawonly) {
  1592.         ads_command(RTSTR, /*MSG0*/"_.LAYER",
  1593.                     RTSTR, /*MSG0*/"_SET", RTSTR, OrbitLayer,
  1594.                     RTSTR, "", RTNONE);
  1595.     }
  1596.  
  1597.     /* Display frame of reference in mode line.  Do this  here  so
  1598.        it's  (a) outside the inner loop, and (b) after changing to
  1599.        the OrbitLayer (if !drawonly) wipes the mode line.  */
  1600.  
  1601.     sprintf(rpt, /*MSG84*/"Reference frame: %s",
  1602.             framei < 0 ? /*MSG85*/"Inertial" :
  1603.              ptab[framei].partname);
  1604.     ads_grtext(-1, rpt, False);
  1605.  
  1606.     /* Initialise last plotted position for each particle. */
  1607.  
  1608.     for (i = 0; i < nptab; i++) {
  1609.         Cpoint(ptab[i].lastpos, ptab[i].position);
  1610.     }
  1611.  
  1612.     while ((nptab > 1) &&
  1613.            (timelimit ? (simtime < endtime) : (n-- > 0)) &&
  1614.            !ads_usrbrk()) {
  1615.         ads_real mindist = 1E100,
  1616.                  maxvel = -1;
  1617.  
  1618.         /* Update acceleration for all bodies. */
  1619.  
  1620.         for (i = 0; i < nptab; i++) {
  1621.             Spoint(ptab[i].acceleration, 0, 0, 0);
  1622.             for (j = 0; j < nptab; j++) {
  1623.                 if (i != j) {
  1624.                     ads_real vdist = ads_distance(ptab[i].position,
  1625.                                                   ptab[j].position),
  1626.                              /* F = G (m1 m2) / r^2 */
  1627.                              force = -GRAVCON *
  1628.                               ((ptab[i].mass * ptab[j].mass) /
  1629.                                (vdist * vdist)),
  1630.                              /* F = ma, hence a = F/m */
  1631.                              accel = force / ptab[i].mass;
  1632.  
  1633.                     mindist = min(mindist, vdist);
  1634.                     if (vdist == 0.0) {
  1635.                         ads_printf(/*MSG86*/"Collision between %s and %s\n",
  1636.                                    ptab[i].partname, ptab[j].partname);
  1637.                     } else {
  1638.                         /* Update vector components of acceleration. */
  1639.                         for (k = X; k <= Z; k++) {
  1640.                             ptab[i].acceleration[k] +=
  1641.                                accel * (ptab[i].position[k] -
  1642.                                         ptab[j].position[k]) / vdist;
  1643.                         }
  1644.                     }
  1645.                 }
  1646.             }
  1647.         }
  1648.  
  1649.         /* Update velocity for all bodies. */
  1650.  
  1651.         for (i = 0; i < nptab; i++) {
  1652.             ads_point pacc;
  1653.             ads_real pvel;
  1654.  
  1655.             addvec(pacc, ptab[i].velocity, ptab[i].acceleration);
  1656.             pvel = sqabsv(pacc);
  1657.             maxvel = max(maxvel, pvel);
  1658.         }
  1659.         maxvel = sqrt(maxvel);
  1660.  
  1661.         deltat = stepsize * (mindist / maxvel);
  1662.         deltat = max(stepmin, deltat);
  1663.         for (i = 0; i < nptab; i++) {
  1664.             sumvec(ptab[i].velocity, ptab[i].velocity,
  1665.                    deltat, ptab[i].acceleration);
  1666.         }
  1667.  
  1668.         /* Update position for all bodies. */
  1669.  
  1670.         for (i = 0; i < nptab; i++) {
  1671.             sumvec(ptab[i].position, ptab[i].position,
  1672.                    deltat, ptab[i].velocity);
  1673.         }
  1674.  
  1675.         /* If we're viewing from one particle's frame of reference,
  1676.            translate the view to adjust  for  our  home  particle's
  1677.            most recent motion.  */
  1678.  
  1679.         if (framei >= 0) {
  1680.             ads_point delta;
  1681.  
  1682.             subvec(delta, ptab[framei].lastpos, ptab[framei].position);
  1683.             for (i = 0; i < nptab; i++) {
  1684.                 addvec(ptab[i].position, ptab[i].position, delta);
  1685.             }
  1686.         }
  1687.  
  1688.         /* Display motion since last update. */
  1689.  
  1690.         for (i = 0; i < nptab; i++) {
  1691.             if (drawonly) {
  1692.                 ads_grdraw(ptab[i].lastpos, ptab[i].position,
  1693.                            ptab[i].colour, False);
  1694.             } else {
  1695.                 if (lastcolour != ptab[i].colour) {
  1696.                     ads_command(RTSTR, /*MSG0*/"_.COLOUR",
  1697.                                 RTSHORT, ptab[i].colour,
  1698.                                 RTNONE);
  1699.                     lastcolour = ptab[i].colour;
  1700.                 }
  1701.                 ads_command(RTSTR, /*MSG0*/"_.LINE", RTPOINT, ptab[i].lastpos,
  1702.                             RTPOINT, ptab[i].position, RTSTR, "", RTNONE);
  1703.             }
  1704.             Cpoint(ptab[i].lastpos, ptab[i].position);
  1705.         }
  1706.  
  1707.         /* Update the step number, simulated time, and display
  1708.            the values on the status line if selected. */
  1709.  
  1710.         stepno++;
  1711.         simtime += deltat;
  1712.         rpt[0] = EOS;
  1713.         if (showtime) {
  1714.             sprintf(rpt, /*MSG87*/"Year %.2f", simtime);
  1715.         }
  1716.         if (showstep) {
  1717.             if (rpt[0] != EOS)
  1718.                 strcat(rpt, "  ");
  1719.             sprintf(rpt + strlen(rpt), /*MSG88*/"Step %ld", stepno);
  1720.         }
  1721.         if (rpt[0] != EOS)
  1722.             ads_grtext(-2, rpt, False);
  1723.     }
  1724.  
  1725.     /* Restore original layer, colour, and UCS. */
  1726.  
  1727.     if (!drawonly) {
  1728.         ads_command(RTSTR, /*MSG0*/"_.LAYER", RTSTR, /*MSG0*/"_SET",
  1729.                     RTSTR, clayer.resval.rstring, RTSTR, "", RTNONE);
  1730.         ads_command(RTSTR, /*MSG0*/"_.COLOUR",
  1731.                     RTSTR, ocolour.resval.rstring, RTNONE);
  1732.     }
  1733.     free(clayer.resval.rstring);
  1734.     free(ocolour.resval.rstring);
  1735.     ads_command(RTSTR, /*MSG0*/"_.UCS", RTSTR, /*MSG0*/"_Prev", RTNONE);
  1736.     CommandE();
  1737.  
  1738.     /* If called as a function, return the simulation end time.
  1739.        If called as a command, print the end time and step count. */
  1740.  
  1741.     if (functional)
  1742.         ads_retreal(simtime);
  1743.     else
  1744.         ads_printf(/*MSG89*/"\n\
  1745. End at step %ld, %.2f years.\n", stepno, simtime);
  1746. }
  1747.  
  1748. /*  SETGRAV  --  Set modal variables.  */
  1749.  
  1750. static void setgrav()
  1751. {
  1752.     struct resbuf rbb;
  1753.     ads_name vname;
  1754.     long l;
  1755.  
  1756.     /* Build the SSGET entity buffer chain to filter for block
  1757.        insertions of the named block on the named layer. */
  1758.  
  1759.     rbb.restype = 2;                  /* Block name */
  1760.     rbb.resval.rstring = ModalBlock;
  1761.     rbb.rbnext = NULL;
  1762.  
  1763.     if (ads_ssget(/*MSG0*/"_X", NULL, NULL, &rbb, vname) == RTNORM) {
  1764.  
  1765.         /* Delete all definitions found. */
  1766.  
  1767.         if (ads_sslength(vname, &l) < 0)
  1768.             l = 0;
  1769.  
  1770.         if (l > 0) {
  1771.             ads_name ename;
  1772.  
  1773.             if (ads_ssname(vname, 0L, ename) == RTNORM) {
  1774.                 CommandB();
  1775.                 ads_command(RTSTR, /*MSG0*/"_.DDATTE", RTENAME, ename, RTNONE);
  1776.                 CommandE();
  1777.                 varset();
  1778.             }
  1779.         }
  1780.         ads_ssfree(vname);            /* Release selection set */
  1781.     }
  1782. }
  1783.  
  1784. /*  UPDATE  --  Update masses to position them at the end of the current
  1785.                 point in the simulation.  */
  1786.  
  1787. static void update()
  1788. {
  1789.     int i;
  1790.  
  1791.     reset();                          /* Reset the simulation */
  1792.     for (i = 0; i < nptab; i++) {
  1793.         ads_name ename;
  1794.  
  1795.         if (ads_handent(ptab[i].partblock, ename) == RTNORM) {
  1796.             struct resbuf *rent, *rb, *vb;
  1797.  
  1798.             memcpy(&pt, &(ptab[i]), sizeof(particle));
  1799.             rent = ads_entget(ename);
  1800.             if (rent != NULL && ((vb = entitem(rent, 10)) != NULL)) {
  1801.                 Cpoint(vb->resval.rpoint, ptab[i].position);
  1802.                 if (ads_entmod(rent) != RTNORM) {
  1803.                     ads_printf(/*MSG90*/"\
  1804. UPDATE: Unable to update position.\n");
  1805.                 }
  1806.             }
  1807.             ads_relrb(rent);
  1808.  
  1809.             /* Update attributes */
  1810.  
  1811.             while (True) {
  1812.                 ads_entnext(ename, ename);
  1813.                 rent = ads_entget(ename);
  1814.                 if (rent == NULL) {
  1815.                     ads_printf(/*MSG91*/"UPDATE: Can't read attribute.\n");
  1816.                     break;
  1817.                 }
  1818.                 rb = entitem(rent, 0);  /* Entity type */
  1819.                 if (rb != NULL) {
  1820.                     if (strcmp(rb->resval.rstring, /*MSG0*/"ATTRIB") != 0)
  1821.                         break;
  1822.                     rb = entitem(rent, 2);  /* Attribute tag */
  1823.                     vb = entitem(rent, 1);  /* Attribute value */
  1824.                     if (rb != NULL && vb != NULL) {
  1825.                         int j;
  1826.  
  1827.                         for (j = 0; j < ELEMENTS(partatt); j++) {
  1828.                             if (strcmp(rb->resval.rstring,
  1829.                                        partatt[j].attname) == 0) {
  1830.                                 /* All right, we've found an attribute in
  1831.                                    the  table.   Edit  the new value into
  1832.                                    it, update the entity,  and  continue. */
  1833.                                 char *sp = vb->resval.rstring;
  1834.                                 char ebuf[32];
  1835.  
  1836.                                 sprintf(ebuf, "%.12g", *partatt[j].attvar);
  1837.                                 vb->resval.rstring = ebuf;
  1838.                                 if (ads_entmod(rent) != RTNORM) {
  1839.                                     ads_printf(
  1840.                                        /*MSG92*/"\
  1841. UPDATE:  Cannot update %s attribute.\n",
  1842.                                        partatt[j].attname);
  1843.                                 }
  1844.                                 /* Restore string buffer */
  1845.                                 vb->resval.rstring = sp;
  1846.                                 break;
  1847.                             }
  1848.                         }
  1849.                     }
  1850.                 }
  1851.                 ads_relrb(rent);
  1852.             }
  1853.         } else {
  1854.             ads_printf(/*MSG93*/"\nCannot retrieve %s.\n", ptab[i].partname);
  1855.         }
  1856.     }
  1857. }
  1858.  
  1859. #ifdef  HIGHC
  1860.  
  1861. /*  Earlier versions of High C put abort() in the same module with exit();
  1862.     ADS defines its own exit(), so we have to define our own abort(): */
  1863.  
  1864. static void abort(void)
  1865. {
  1866.     ads_abort("");
  1867. }
  1868.  
  1869. #endif  /* HIGHC */
  1870.